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INTRODUCTION. 

The  horseshoe-vortex  and  associated  corner  flows  which  occur  when  a  hlunt 
obstruction  is  placed  within  an  approaching  boundary  layer  represent  a 
fundamental  three-dimensional  viscous  flow  of  considerable  interest  and 
importance.  Examples  of  this  type  of  flow  include  the  flows  near  an  aircraft 
wing/fuselage  junction  and  near  a  submarine  hull/sail  junction.  The  feature 
common  to  all  horseshoe  vortex  flows  is  that  a  non-uniform  velocity  in  an 
approaching  boundary  layer  meets  a  local  region  of  adverse  pressure  gradient 
due  to  the  blockage  effect  of  the  obstruction.  This  causes  a  three- 
dimensional  boundary  layer  separation  and  the  formation  of  one  or  more 
horseshoe  vortices  around  the  obstruction.  In  addition  to  the  leading  edge 
region,  the  associated  corner  flows  downstream  of  the  leading  edge  are  also  of 
interest,  since  they  contain  streamwise  vortices  which  affect  the  performance 
of  both  the  airfoi-  or  strut  and  also  other  devices  located  downstream. 

The  problem  of  horseshoe-vortex/corner  flow  has  been  investigated 
previously  [1-2]  by  numerical  solution  of  the  Navier-Stokes  equations  using  a 
linearized  block  implicit  (LBI)  scheme  with  Douglas-Gunn  splitting.  Solutions 
for  laminar  flow  past  an  elliptical  leading  edge  mounted  normal  to  a  flat  plate 
endwall  have  been  computed  for  both  zero  and  five-degree  angles  of  incidence, 
with  chordal  Reynolds  number  of  400,  and  Mach  number  of  0.2.  Turbulent  flow 
cases  have  been  computed  for  both  unswept  and  45-degree  swept  elliptical 
leading  edges  mounted  on  a  flat  plate,  with  Reynolds  number  of  310,000  and  Mach 
number  of  0.05.  More  recently,  Briley,  Buggeln  and  McDonald  [3]  have  found 
that  appropriate  inf low/outflow  boundary  conditions  provide  conditional 
stability  and  rapid  convergence  for  the  three-dimensional  algorithm  and  have 
considered  alternative  forms  of  artificial  dissipation.  An  accurate  (grid 
Insensitive)  solution  was  computed  for  a  laminar  horseshoe  vortex  flow  using 
only  small  amounts  of  artificial  dissipation.  Very  good  convergence  rates  were 

_ p 

obtained  for  this  laminar  flow  case,  with  residuals  reduced  by  a  factor  of  10 
in  30  and  60  iterations,  respectively,  for  15x15x15  and  30x30x30  grids.  This 
work  is  extended  In  the  present  study  to  address  turbulent  flows  using  both 
two-equation  (k-e )  and  algebraic  mixing  length  turbulence  models,  with  viscous 
sublayer  resolution. 

Recent  experimental  work  on  turbulent  horseshoe  vortex  flows  includes  that 
of  McMahon,  Hubbart  and  Kubendran  [4],  Moore  and  Forlini  [5],  Eckerle  and 
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Langston  (6]  and  Dickinson  [7].  Moore  and  Moore  [8)  have  performed  numerical 
computations  using  the  inviscid  Ruler  equations  for  a  horseshoe  vortex  flow  and 
have  compared  these  results  with  the  measurements  of  Moore  and  Forlini  (5). 
Calculations  for  a  hlunt-fin  induced  shock  wave  and  boundary  layer  Interaction 
flow  containing  a  horseshoe  vortex  have  been  reported  by  Hung  and  Kordulla  [9]. 

In  the  present  report,  the  methods  used  in  solving  the  three-dimensional 
Navier-Stokes  equations  and  two-equation  turbulence  model  are  first  reviewed. 
Turbulence  models  based  on  solution  of  the  k-e  equations  and  on  an  algebraic 
mixing  length  model  are  then  implemented  within  a  vectorized  code  using 
coordinate  distributions  which  provide  resolution  of  viscous  sublayer  regions. 
Resolution  of  the  viscous  sulayer  region  is  believed  to  be  important  for 
three-dimensional  flows,  for  which  wall-function  boundary  conditions  are  of 
limited  value.  Test  calculations  are  then  performed  for  a  two-dimensional 
turbulent  flow  to  develop  the  computational  methodology  necessary  to  obtain 
rapid  convergence.  Finally,  three-dimensional  turbulent  horseshoe  vortex  flows 
are  computed  using  both  k-s  and  mixing  length  models.  These  results  are 
displayed  graphically  and  are  evaluated  by  comparison  with  recent  detailed 
experimental  measurements  obtained  by  McMahon,  Hubbartt  and  Kubendran  [4]. 

GOVERNING  EQUATIONS 

The  three-rtimensional  compressible  Navier-Stokes  equations  are  solved  here 
for  low  Mach  number  and  with  an  assumption  of  constant  stagnation  enthalpy. 

For  these  conditions,  steady  flow  solutions  closely  approximate  an 
incompressible  constant  density  flow  (cf.  Briley,  McDonald  and  Shamroth  [10]). 

A  zonal  approach  is  used  wherein  the  flow  is  computed  only  in  a  subregion  of 
the  overall  flow  field,  near  the  leading  edge  (Fig.  1).  The  form  of  the 
governing  equations  solved  permits  the  use  of  general  nonorthogonal  body-fitted 
coordinate  systems,  and  is  obtained  by  a  transformation  from  Cartesian  to 
general  nonorthogonal  coordinates.  The  Cartesian  velocity  components  uj  and 
density  p  are  retained  as  dependent  variables  in  the  transformed  system  of 
equations.  The  pressure,  temperature  and  stagnation  enthalpy  are  denoted  p,  T 
and  hn,  respectively. 

All  variables  are  nondimens  tonal  in  the  present  formulation,  having  been 
normalized  by  reference  quantities  denoted  by  a  subscript  'r'.  The  quantities 
Pr,  Ur ,  Tr  and  Lr  denote  reference  values  for  density,  velocity,  temperature 
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and  length,  respectively.  The  reference  pressure,  enthalpy  and  time  are  taken 
as  prlJr,  CpTr  and  Lr/Ur,  respecti vely,  where  cp  Is  the  specific  heat  at 
constant  pressure.  The  specific  heat  ratio  is  y»  end  Mr  ■  Ur/cr  is  a  reference 
Mach  number,  where  cr  is  the  reference  sound  speed  defined  by  cr  =  yRTr, 
and  R  is  the  gas  constant.  The  reference  Reynolds  number  Re  is  defined  by 
prUrLr/ur,  where  ur  is  a  reference  viscosity. 

The  transformation,  T,  from  Cartesian  coordinates  xj  to  computational 
coordinates  yd  is  given  by 

T  =  yJ  (xp  i,j  =  1,2,3  (1) 

Spatial  derivatives  are  transformed  according  to 


_3 

3x 


(2) 


where,  unless  otherwise  stated,  the  summation  convention  is  used  for  repeated 

indices,  and  y1 .  =  3yJ/3x  .  The  coordinate  system  is  defined  by  specifying 

,1  1 

the  Cartesian  coordinates  of  each  computational  grid  point.  The  partial 
derivatives  3x^/3)^  of  the  inverse  transformation  T  *  =  x^ (y^ )  are  then 
computed  using  three-point  second-order  difference  formulas  with  uniform 
spacing  of  the  computational  coordinates  yd.  For  convenience,  the  yJ 
coordinates  are  normalized  to  give  a  unit  mesh  spacing  Ayd  =  1  for  each 
coordinate.  The  transformation  derivatives  3yd / 3xj  are  then  computed  from 
3xj/3yd  using  standard  procedures  for  computing  derivatives  of  inverse 
functions  (cf.  Kaplan  [11]). 

The  transformed  Navier-Stokes  equations  can  be  written  in  the  following 
nondimensional  form:  The  continuity  equation  is 


3p 

3t 
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The  kth  component  of  the  momentum  equation  is  given  by 
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where  <5^  is  the-  Kronecker  delta  function.  The  shear  stress  is  given 

by 
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The  equation  of  state  and  definition  of  stagnation  enthalpy  can  be  expressed 
for  a  perfect  gas  as 

P  -  pt/ymJ  (6) 

hQ  =  T/(Y-l)Mr2  +  q2/2  (7) 

where  q2  =  S^u^uj.  Although  it  is  not  necessary,  it  is  both  convenient 
and  computationally  worthwhile  for  the  present  problem  to  assume  that  h0  is  a 
constant  and  to  omit  solution  of  the  energy  equation.  This  results  in 
negligible  error  for  steady  flow  at  low  Mach  number  with  no  heat  addition. 
Equations  (6)  and  (7)  can  then  be  combined  to  produce  an  adiabatic  equation  of 
state 


p  =  p (h  -  q2/2)  (y-1)/y 

°  (8) 

which  is  used  to  eliminate  pressure  as  a  dependent  variable  in  Eq.  (4). 

Two  turbulence  models  have  been  utilized  in  the  present  investigation. 
These  are  a  two-equation  (k-e)  model  based  on  that  due  to  Jones  and  Launder 
[121,  and  a  simple  algebraic  mixing  length  model.  Each  of  these  models  is 
formulated  to  provide  resolution  of  the  viscous  sublayer  region  near  solid 
walls.  The  turbulent  viscosity  is  assumed  to  be  isotropic  and  is  added  to  the 
laminar  viscosity  to  obtain  a  total  effective  viscosity,  pe  =  u  +  p-p. 

The  governing  turbulence  kinetic  energy  (k)  and  dissipation  rate  (e)  are 
given  by  Launder  and  Spalding  [13].  Applying  the  general  coordinate 
transformation,  these  equations  can  be  written  as 
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Jones  and  Launder  (12]  have  suggested  adding  additional  terms  to  the  turbulence 
energy  and  dissipation  equations  for  use  in  computing  low  Reynolds  number 
turbulent  flows.  An  extra  term  was  added  to  the  k  equation  for  computational 
reasons  to  counteract  the  effect  of  applying  e  =  0  as  a  wall  boundary 
condition.  In  the  present  study,  these  extra  terms  were  omitted  and  a 
zero-gradient  condition  for  e  was  applied  as  a  wall  condition. 

An  algebraic  mixing  length  model  was  used  both  for  flow  predictions  and  as 
a  method  of  initializing  the  flow  field  and  starting  the  k-c  model.  In  the 
mixing  length  model,  the  turbulent  viscosity  is  obtained  from 


,  1/2 
UT  =  pH 2  (20^0^) 


It  remains  to  specify  a  length  scale  distribution  appropriate  for  the  problem 
under  consideration.  The  horseshoe-vortex/corner  flow  of  interest  here  has 
moderately  thin  shear  layers  on  the  endwall  and  strut  surface,  and  the  length 
scale  distribution  is  thus  adapted  from  previous  turbulence  models  for 
turhulent  boundary  layers  and  taken  to  be  the  conventionally  defined  mixing 
length  of  Prandtl.  The  distribution  of  mixing  length  given  by  McDonald  and 
Fish  (14]  has  proven  effective  for  a  wide  range  of  two-dimensional  turbulent 
boundary  layers  and  is  easily  adapted  for  present  use.  This  distribution  is 
given  by 

£  =  ^  £oo  tanh  (<d/£oo)  (17) 

where  £  is  mixing  length,  £co  is  an  outer-region  value  of  mixing  length,  d  is 
distance  from  the  wall,  <  is  the  von  Karman  constant  (taken  here  as  0.4),  and 
is  the  Van  Driest  sublayer  damping  function  given  by 


=  1  -  exp  (-y  /A  ) 


where 
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and  u  =  (t  /p)  and  T  is  the  wall  shear  stress.  The  constant  A  is  taken  as 
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26.0.  For  equilibrium  boundary  layers,  the  freestream  mixing  length,  is 
observed  to  have  a  constant  value  of  about  0.1  6,  where  6  is  the  local  boundary 
layer  thickness. 

The  length  scale  distribution  of  Eq.  (17)  is  adapted  for  present  use  by 
taking  d  as  distance  to  the  nearest  wall  and  by  assigning  la,  a  distribution 
based  on  two-dimensional  momentum  integral  estimates  of  the  boundary  layer 
growth  rates  on  the  endwall.  A  constant  value  of  freestream  mixing  length  was 
used  for  the  strut  surface.  Two  boundary  layer  thickness  distributions,  for 
the  endwall  boundary  layer  and  62  for  the  strut  surface,  are  obtained  and  an 
outer  or  maximum  mixing  length  scale  associated  with  each  boundary  layer  is 
determined  from  the  formula  la,  =  0.1  6.  In  the  core  region  outside  the 
boundary  layer  the  length  scale  is  exponentially  damped  to  a  negligible  value 
over  a  distance  of  about  two  boundary  layer  thicknesses.  In  the  overlap  corner 
region,  the  outer  or  maximum  length  scale,  la,,  is  taken  as  the  greater  of  the 
two  values  associated  with  the  strut  and  endwall  boundary  layers. 


METHOD  OF  SOLUTION 

The  basic  algorithm  considered  here  has  been  described  by  Briley  and 
McDonald  [IS,  16]  and  employs  a  formal  time  linearization  to  produce  a 
noniterative  fu 1 Ly-coupled  approximation  for  nonlinear  systems  of  equations, 
which  is  solved  in  block-implicit  form  using  an  ADI  scheme  with  consistent 
intermediate  steps.  For  a  linear  scalar  diffusion  equation,  this  algorithm 
reduces  to  a  classical  ADI  scheme  considered  by  Douglas  and  Gunn  fl7).  Warming 
and  Beam  [18,  19]  have  introduced  a  very  concise  derivation  of  this  same 
algorithm  using  approximate  factorization  of  the  linearized  approximation 
written  in  'delta'  form.  The  works  of  Pulliam  and  Steger  [20],  Thomas  and 
Lombard  [21]  and  Shamroth,  McDonald  and  Briley  [22]  are  representative  of 
numerous  investigations  which  have  employed  this  basic  algorithm. 

Linearization  and  Time  Differencing 

The  nonlinear  system  of  governing  equations  is  first  written  (at  a  single 
grid  point)  in  the  following  form: 


3H(4>)/3t  =  D(4>)  +  S  (4>) 


where  <t>  is  the  column-vector  of  dependent  variables,  H  and  S  are  column-vector 
algebraic  functions  of  4> »  and  D  is  a  column  vector  whose  elements  are  the 
spatial  differential  operators  which  generate  all  spatial  derivatives  appearing 
in  the  governing  equation  associated  with  that  element. 

The  solution  procedure  is  based  on  the  following  two-level  implicit 
time-difference  approximation  of  (20): 


(Hn+1  -  Hn)/At  =  B(Dn+1  +  Sn+1)  +  (1-6)  (Dn  +  Sn) 


where,  for  example,  Hn+  denotes  H(<J>n+  )  and  At  =  tn+l  -  tn.  The  parameter 
6  (0.5  <  B  <  1)  permits  a  variable  time-centering  of  the  scheme,  with  a 
truncation  error  of  order  f  At  ,  (8  -  1/2  At)]. 


A  local  time  linearization  (Taylor  expansion  about  $n)  of  requisite 
formal  accuracy  is  introduced,  and  this  serves  to  define  a  linear  differential 
operator,  L,  such  that 

Dn+l  =  D°  +  Ln  (0n+1  -  4>n)  +  0  (At2)  (22a) 

Similarly, 

Hn+1  =  Hn  +  (3H/3<»)n  (4>n+1  -  4.")  +  0  (At2)  (22b) 

sn+1  =  sn  +  (as/ 3<t>) n  (4»n+1  -  $n)  +  o  (At2)  (22c) 

Equations  (22a-c)  are  inserted  into  Eq.  (21)  to  obtain  the  following  system 
which  is  linear  in  (f>n+l 


(A  -  BAt  Ln)  (4>n+1  -  cj,n)  =  At  (Dn  +  Sn)  (23) 

and  which  is  termed  a  Linearized  block  implicit  (LBI)  scheme.  Here,  A  denotes 
a  square  matrix  defined  by 

A  H  OH/a<t>)n  -  BAt  Os/3(J))n  (24) 

Equation  (22)  has  0  (At)  accuracy  unless  H  =  4>,  in  which  case  the  accuracy  is 
the  same  as  Eq .  (21). 

Special  Treatment  of  Diffusive  Terms 


Spatial  cross-derivatives  are  present  in  viscous  terms  and  in  added 
artificial  dissipation  terms  of  the  present  formulation,  and  these  cross 
derivative  terms  are  evaluated  explicitly  at  tn.  To  preserve  notational 
simplicity,  it  is  understood  that  all  cross-derivative  terms  appearing  in  Ln 
are  neglected  but  are  retained  in  Dn.  In  addition,  although  diffusion 
coefficients  in  viscous  and  dissipation  terms  are  generally  functions  of  the 
dependent  variables,  these  coefficients  are  not  linearized  and  instead  are 
evaluated  implicitlv  at  tn  during  each  time  step.  Notat i ona  1  ly ,  this  is 
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equivalent  to  neglecting  derivatives  of  these  coefficients  with  respect  to  $  In 
Ln,  which  are  formally  present  In  the  Taylor  expansion  (22a),  but  otherwise 
retaining  all  terms  in  both  Ln  and  Dn. 

It  Is  Important  to  note  that  neglecting  terms  In  Ln  has  no  effect  on 
steady  solutions  of  Eq.  (23),  since  <j>n+l  —  (J>n  =  0  and  thus  Eq.  (23)  reduces 
to  the  steady  form  of  the  equations:  Dn  +  Sn  =  0.  Aside  from  stability 
considerations,  the  only  effect  of  neglecting  terms  in  Ln  is  to  introduce  an 
0(At)  truncation  error. 

Consistent  Splitting  of  the  LBI  Scheme 

To  obtain  an  efficient  algorithm,  the  linearized  system  (23)  is  split 
using  ADI  techniques.  To  obtain  the  split  scheme,  the  multidimensional 
operator,  L,  is  rewritten  as  the  sum  of  three  "one-dimensional”  sub-operators 
Lj_  (i  =  1,  2,  3)  each  of  which  contains  all  terms  having  derivatives  with 
respect  to  the  i-th  spatial  coordinate.  The  split  form  of  Eq.  (23)  can  be 
derived  either  by  following  the  procedure  described  by  Douglas  and  Gunn  in 
their  generalization  and  unification  of  scalar  ADI  schemes,  or  using 
approximate  factorization.  In  either  case,  for  the  present  system  of  equations 
the  split  algorithm  is  given  by 

(A  -  BAth")  (4>*  -  4>°)  =  At  (Dn  +  Sn)  (25a) 

(A  -  BAtL")  ($**  -  4>n)  =  A  ($*  -  4>n)  (25b) 

(A  -  BAtL”)  (<J>n+1  -  4>n)  =  A  (4>**  -  4>n)  (25c) 

where  <p*  and  <)>**  are  consistent  intermediate  solutions.  If  spatial  derivatives 
appearing  in  and  D  are  replaced  by  three-point  difference  formulas,  then 
each  step  in  Eqs.  (25a-c)  can  be  solved  by  a  block-tridiagonal  elimination. 

Combining  Eqs.  (25a-c)  gives 


(A  -  BAtLj)  A  1  (A  -  B'AtL")  A  1  (A  -  BAtL")  (<j>n+1  -  $n) 


=  At  (Dn  +  Sn) 


-I  1- 


(26) 


which  approximates  the  unsplit  scheme  (22)  to  0  (At  ).  Since  the  intermediate 
steps  are  also  consistent  approximations  for  Eq.  (22),  physical  boundary 
conditions  can  be  used  for  and  $**.  Finally,  since  the  Lj  are  homogeneous 
operators,  it  follows  from  Eq'  .  (25a-c)  that  steady  solutions  have  the  property 
that  <f>n  +  l  =  <J>*  =  $**  =  <j>n  and  satisfy 


Dn  +  Sn 


=  0 


(27) 


The  steady  solution  thus  depends  only  on  the  spatial  difference  approximation 
used  for  (27),  and  does  not  depend  on  the  solution  algorithm  itself. 

In  applying  the  split  I.BI  scheme  to  the  turbulence  kinetic  energy 
equation,  it  should  be  noted  that  the  linearization  employed  by  the  viscous 
terms  (turbulent  viscosity  is  evaluated  explicitly  at  tn)  decouples  the 
turbulence  equations  from  the  continuity  and  momentum  equations.  This  means 
that  the  two  turbulence-model  equations  can  be  solved  as  a  separate  2x2  coupled 
system,  once  the  continuity  and  momentum  equations  have  been  solved.  The 
turbulence  kinetic  energy  and  dissipation  equations  contain  highly  nonlinear 
sc  :ce  terms  (RHS  of  Eqs.  9-10)  which  require  linearization  for  implementation 
of  the  LBT  scheme.  These  source  terms  are  first  rewritten  utilizing  Eqs.  (11) 
and  (13)  to  give 


SU)  = 


2 Re  1UTDijDij  ~  Pe 
CuC^jD^  pk  -  C2PE2/k 


(28) 


The  last  term  in  Equation  (11)  is  negligible  for  low  Mach  number  flows  and  has 
been  neglected.  These  source  terms  are  then  linearized  as  follows: 


0 


-P 


as 

3<j) 


CuC iD| j  D( ^  p  + 


-2C2p(i) 

k 


(29) 
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COMPUTED  RESULTS  FOR  A  TURBULENT  HORSESHOE  VORTEX  FLOW 


Flow  Conditions 


Solutions  are  presented  here  for  turbulent  flow  at  zero  Incidence  past  an 
elliptical  leading  edge  geometry  mounted  between  parallel  flat  plate  endwalls. 
The  leading-edge  geometry  for  the  present  calculations  is  the  same  as  that 
considered  by  McMahon,  Hubbartt  and  Kubendran  [4]  for  turbulent  flow 
conditions.  The  flow  geometry  (Fig.  1)  consists  of  a  strut  of  constant 
thickness,  W,  having  an  elliptical  leading  edge  with  1.5:1  ratio  of  major  to 
minor  axis.  The  strut  is  mounted  normal  to  parallel  flat  plate  endwalls  whose 
separation  distance  is  5.0W,  and  whose  leading  edges  are  located  a  distance 
6.0W  upstream  of  the  leading  edge  of  the  strut.  The  length  L  of  the  strut 
within  the  computational  domain  is  2.5W.  The  flow  considered  has  a  Reynolds 
number  Re  =  57,000  (based  on  L)  and  Mach  number  Mr  =  0.05,  each  based  on 
upstream  flow  conditions. 

Boundary  Conditions 

Since  the  computational  domain  is  chosen  to  be  a  region  in  the  immediate 
vicinity  of  the  leading-edge/corner  flow  geometry  (cf.  Fig.  1)  embedded  within 
a  larger  overall  flow  system,  inflow  and  outflow  boundary  conditions  which 
adequately  model  the  interface  between  the  computed  flow  and  the  remainder  of 
the  flow  svstem  are  required.  The  inflow/outflow  conditions  used  are  derived 
from  an  assumed  flow  structure  and  are  chosen  to  provide  inflow  with  prescribed 
stagnation  pressure  (and  stagnation  enthalpy)  in  an  inviscid  core  region  and 
with  a  given  axial  velocity  profile  shape  in  the  endwall  boundary  layer,  and  to 
provide  outflow  with  a  prescribed  distribution  of  static  pressure  in  the  cross 
section.  This  approach  to  inflow/outflow  boundary  conditions  has  been 
discussed  previously  in  [1-2]  and  is  only  summarized  here. 

First,  the  houndary  laver  thickness  6(xj)  on  the  endwall  flat  plate  is 
approximated  by  its  distribution  from  the  Blasius  flat  plate  solution.  At  the 
inflow  boundary,  a  "two-laver"  boundary  condition  is  employed  such  that 

stagnation  pressure,  p0 ,  is  fixed  at  the  free  stream  reference  value  in  the 

3  3 

core  flow  region  (y  >  6)  and  an  axial  velocity  profile  shape  »j/ue  =  f ( y  / 6 ) 
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is  fixed  within  the  boundary  layer  region  (y  <  <S).  Here,  ue  is  the  local 
edge  velocity  which  varies  with  time  and  is  adjusted  after  each  time  step  to 
the  value  consistent  with  pQ  and  the  local  edge  static  pressure,  which  is 
determined  as  part  of  the  solution.  The  remaining  inflow  conditions  are 
u2  =  3  u3/3n  =  3  p/3n  =  0,  where  n  denotes  the  normal  computational 

coordinate,  y1.  For  outflow  conditions,  a  constant  static  pressure  is  imposed 
and  second  derivatives  of  each  velocity  component  are  set  to  zero.  At  no-slip 
surfaces,  each  velocity  component  u^  is  set  to  zero,  and  the  remaining 
condition  applied  to  these  surfaces  is  that  the  derivative  of  pressure  in  the 
direction  normal  to  the  surface  is  zero.  This  condition  approximates  the 
normal  momentum  equation  to  order  Re-1  for  viscous  flow  at  a  no-slip  surface. 
The  final  boundary  to  be  considered  is  the  plane  parallel  to  the  endwall  and  in 
the  free  stream.  This  boundary  is  assumed  to  be  a  plane  of  symmetry,  so  that 
the  flow  represented  is  that  past  the  strut  mounted  between  parallel  flat 
plates. 

Although  a  uniform  grid  was  used  for  the  transformed  coordinates  y j ,  the 
grid  used  is  highly  nonuniform  in  physical  coordinates  x^,  as  shown  in 
Fig.  1,  and  was  chosen  to  provide  resolution  of  several  length  scales  known  to 
be  present  for  this  type  of  flow.  Care  was  taken  to  provide  resolution  of  the 
boundary  layers  and  viscous  sublayers  on  the  strut  and  endwall  surfaces  of  the 
shear  layer  near  the  leading-edge  stagnation  line,  and  of  the  corner  flow 
region  very  near  the  endwall/leading-edge  intersection. 

A  pseudo-time  dependent  iteration  procedure  was  employed  in  which  the 
square  matrix  A  in  the  LBI  scheme  (25)  is  replaced  by  a  modified  matrix  A  <•  AB 
where  B  is  a  diagonal  conditioning  matrix.  This  technique  was  employed  bv 
Briley,  McDonald  and  Shamroth  [10]  to  improve  the  convergence  rate  of  this 
algorithm  when  applied  to  low  Mach  number  flows.  In  the  present  calculations, 
a  diagonal  matrix  B  whose  diagonal  elements  are  (l/yMr2,  1,  1,  1)  was  used. 

In  addition,  a  spatially  varying  time  step  was  used  to  avoid  instability 
analogous  to  that  encountered  for  the  three-dimensional  scalar  convection 
equation,  and  to  improve  the  convergence  rate. 


Computed  Results 

Test  calculations  were  performed  for  a  two-dimensional  turbulent  flow  to 
develop  the  computational  methodology  necessary  to  obtain  rapid  convergence 
using  the  two-equation  turbulence  model  and  coordinate  distributions  which 


provide  resolution  of  viscous  sublayer  regions.  The  test  problem  has  the  same 
geometry  and  flow  conditions  as  the  horseshoe  vortex  flow  to  be  considered 
subsequently,  but  is  restricted  to  the  two-dimensional  region  away  from  the 
strut  endwali  Intersection  where  the  horseshoe  vortex  forms.  The  flow 
conditions  for  these  test  calculations  are  Reynolds  number  Re  =  57,000  based  on 
body  width  and  inflow  Mach  number  Mr  =  0.05.  A  60x50  grid  was  used,  with  a 
coordinate  stretching  such  that  the  mesh  spacing  adjacent  to  the  wall  is 
approximately  0.0005.  The  coordinates  and  steady  flow  solution  for  these 
calculations  are  shown  in  Fig.  2. 

The  Mach  number  contours  in  Fig.  2  give  an  indication  of  the  thin  boundary 
layer  associated  with  turbulent  flow  at  this  Reynolds  number.  Plots  of  maximum 
normalized  residual  and  maximum  increment  of  solution  variable  versus  iteration 
number  are  shown  in  Figs.  3-4  for  both  mixing  length  and  k-e  turbulence  models. 
A  dissipation  parameter  (a  =  0.1)  was  used  for  each  of  these  calculations.  The 
use  of  larger  values  of  dissipation  did  not  have  a  substantial  effect  on  the 
convergence  behavior.  It  should  be  noted  that  for  scalar  convection  in  one 
dimension,  the  value  a  =  1.0  is  equivalent  to  two-point  upwind  differencing, 
and  thus  the  value  o  =  0.1  is  relatively  small.  In  the  mixing  length 
calculations,  the  residual  was  reduced  by  a  factor  of  about  105  in  100 
iterations.  The  k-e  model  was  started  after  25  iterations  using  the  mixing 
length  model,  and  although  the  convergence  was  not  as  fast  as  using  the  mixing 
length,  the  residual  was  reduced  by  about  three  orders  of  magnitude  in  100 
iterations.  It  is  not  clear  whether  the  difference  in  convergence  behavior 
between  these  two  cases  is  due  to  the  physical/mathematical  properties  of  the 
k-e  turbulence  model  or  to  details  of  the  numerical  solution  procedure.  These 
calculations  required  about  27  seconds  of  CRAY-XMP  CPU  time  using  a  vectorized 
code . 

Solutions  for  three-dimensional  horseshoe  vortex  flow  were  computed  using 
both  mixing  length  and  k-e  turbulence  models.  The  geometry  and  flow  conditions 
correspond  to  those  measured  by  McMahon,  Hubbartt  and  Kubendran  (Ref.  4).  The 
Reynolds  number  is  Re  =  57,000  and  a  near  incompressible  Mach  number  Mr  =  0.05 
was  used.  The  endwali  boundary  layer  has  a  thickness  of  about  0.45  W  at  a 
distance  of  3  W  upstream  of  the  leading  edge.  A  60x50x40  grid  (120,000  points) 
was  used,  with  mesh  spacings  of  0.00075  W  and  0.0005  W  adjacent  to  the  endwali 
and  strut  surfaces,  respectively.  The  overall  computational  domain  extended  to 
a  distance  of  3.5  W  from  the  endwali  and  strut  surfaces,  where  symmetry 


conditions  were  applied.  Allowing  for  symmetry,  the  overall  domain  (7  W  x  7  W) 
Is  smaller  than  the  wind  tunnel  cross  section  (18.5  W  x  18.8  W)  In  Ref.  4,  and 
hence  there  Is  a  stronger  blockage  effect  In  the  computed  results  due  to  the 
body  in  the  region  downstream  of  the  leading  edge.  The  k-e  model  was  started 
after  100  iterations  using  the  mixing  length  model,  and  convergence  to  three 
orders  of  residual  reduction  required  about  200  Iterations  in  three 
dimensions.  In  the  three-dimensional  k-e  calculations  it  was  found  necessary 
to  use  a  larger  dissipation  parameter  (a  =  1.0)  than  in  the  mixing  length 
calculations  (o  =  0.1).  These  solutions  required  about  2.5  hours  of  CRAY-XMP 
CPU  time  on  a  vectorized  code,  using  120,000  grid  points. 

The  computed  flow  behavior  using  the  k-e  turbulence  model  is  first 
discussed  and  compared  with  experimental  measurements.  A  comparison  of  the  k-e 
and  mixing  length  solutions  will  be  given  subsequently.  Some  computer¬ 
generated  flow  visualization  results  are  given  in  Figs.  5a-b.  The  computation¬ 
al  grid  shown  is  the  grid  used  for  the  solution,  except  that  approximately  one- 
half  of  the  grid  points  in  each  direction  have  been  omitted.  The  contours  of 
surface  pressure  show  the  three-dimensional  behavior,  with  low  pressure  near 
the  leading-edge/endwall  intersection.  The  particle  path  lines  shown  are 
useful  as  a  visualization  of  the  three-dimensional  horseshoe  vortex  flow 
structure.  Two  groups  of  particle  lines  are  shown.  The  first  group  includes 
particle  lines  which  originate  at  the  upsteam  inflow  boundary  in  a  line  of 
points  parallel  to  the  endwall  and  located  very  near  the  endwall.  These 
particle  lines  remain  very  close  to  the  endwall  but  are  deflected  at  the 
saddle-point  separation  line  and  remain  outside  the  horseshoe  vortex.  The 
second  group  of  particle  lines  originates  from  a  vertical  line  of  points 
Located  near  the  plane  of  symmetry  and  just  upstream  of  the  leading  edge. 

These  particle  lines  originate  behind  the  saddle-point  separation  and  show  the 
downward  flow  near  the  intersection  of  the  leading  edge  and  endwall  surfaces. 
Some  of  these  particle  lines  near  the  endwall  flow  in  a  direction  opposite  the 
freestream  and  toward  the  saddle-point  separation  line  and  then  remain  neat 
this  separation  line  downsteam  of  the  leading  edge.  The  remaining  particle 
lines  stav  much  closer  to  the  body  surface.  The  overall  flow  behavior  revealed 
bv  these  computer-generated  particle  lines  is  verv  similar  to  that  observed  in 
numerous  flow  visualization  experiments  for  horseshoe  vortex  flows  of  this 
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The  computed  flow  near  the  strut/endwal 1  Intersection  Is  compared  with 
surface  oil-flow  visualization  photographs  In  Fig.  6.  Vector  velocity  plots  In 
the  plane  of  grid  points  adjacent  to  the  no-slip  endwall  surface  are  shown, 
together  with  the  selected  streamlines.  A  saddle-point  flow  separation  and 
reverse  flow  region  are  clearly  visible  In  both  computed  and  experimental 
results  and  there  is  good  agreement  In  the  location  of  this  separation  point. 
The  experimental  flow  visualization  has  a  second  saddle-point  separation 
located  near  the  leading  edge.  The  computed  flow  also  has  a  second 
saddle-point  separation,  although  It  is  located  much  closer  to  the  leading  edge 
than  that  evident  In  the  flow  visualization  and  Is  not  clearly  visible  In  Fig. 


McMahon,  Hubbartt  and  Kubendran  (Ref.  4)  measured  streamwlse  and 
transverse  velocity  components  at  streamwlse  stations  both  upstream  and 
downstream  of  the  strut  leading  edge,  as  Indicated  In  Fig.  7.  The  computed 
mean  flow  velocities  at  three  of  the  measured  stations  are  compared  with  these 
measurements  in  Figs.  8a-c.  The  comparison  at  station  1,  located  11  mm 
upstream  of  the  leading  edge,  Is  shown  in  Fig.  8a.  This  station  Is  located 
within  the  reversed  flow  region  of  the  saddle-poirt  flow  separation  and  a  small 
reversed  flow  region  can  be  seen  in  the  computed  streamwlse  velocity  contours. 
In  general,  there  is  good  agreement  between  computed  and  measured  and  mean  flow 
velocities.  The  computed  and  measured  transverse  velocity  vectors  generally 
agree  as  to  strength,  direction  and  skewing  near  the  wall.  The  computed  and 
measured  velocities  at  stations  4  and  5,  located  in  the  corner  region 
downstream  of  the  leading  edge,  are  shown  In  Figs.  8h  and  8c.  Although  the 
strength  of  the  secondary  flow  is  approximately  correct,  the  details  of  the 
corner  vortex  and  the  distortion  of  the  streamwlse  velocity  contours  by  the 
corner  vortex  are  not  accurately  predicted.  It  should  be  noted  that  the 
computed  flow  has  increased  blockage  effects  since  the  computed  flow  region  Is 
smaller  than  the  wind  tunnel  cross  section,  as  mentioned  earlier.  It  should 
also  be  noted  that  the  corner  vortex  is  relatively  weak,  with  tranverse 
velocities  on  the  order  of  5  percent  of  the  streamwise  velocity.  It  seems 
Likely  that  the  k-e  turbulence  model  used  Is  inadequate  for  the  corner  flow 
region,  although  numerical  truncation  error  may  also  affect  this  comparison 
despite  the  relatively  high  grid  resolution  in  this  calculation. 
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Some  detail*?  of  the  computed  corner  secondary  flow  at  station  4  are  shown 


in  Fig.  9.  Transverse  veLocitv  vector  plots  are  given  for  three  different 
Levels  of  magnification  in  the  corner  region.  The  vector  magnitudes 
a re  renormalized  for  each  plot  and  thus  indicate  only  flow  direction  and 
relative  magnitude  within  each  plot.  It  can  he  seen  from  Fig.  9  that  there  is 
a  triple  vortex  structure  in  the  computed  corner  flow.  The  vortices  are 
counter-rotat ing  and  have  decreasing  strength  as  the  corner  is  approached. 
Typical  velocities  in  the  three  plots  shown  are  of  the  order  of  0.1,  0.01  and 
0.0005  times  the  freestream  velocity.  Finally,  a  comparison  of  computed 
results  obtained  using  both  a  k-e  and  mixing  length  turbulence  models  is  given 
in  Fig.  10.  The  two  solutions  shown  are  identical  in  every  respect  except  for 
the  turbulence  model  used.  The  computed  streamwise  and  transverse  velocity 
fields  are  shown  in  Fig.  10  for  station  4,  which  is  downstream  of  the  leading 
edge.  Although  some  difference  between  the  two  solutions  is  evident  in  the 
streamwise  velocity  contours  very  near  the  corner,  the  transverse  velocity  away 
from  the  corner  is  not  affected  by  this  difference,  ar d  overall  the  two 
turbulence  models  give  similar  results. 

CONCLUDING  REMARKS 

Solutions  of  the  three-dimensional  Reynolds-averaged  Navier-Stokes 
equations  have  been  computed  for  turbulent  flow  past  a  strut/endwall 
configuration  for  which  experimental  measurements  are  available.  Computed 
results  were  obtained  for  both  two-equation  (k-e)  and  mixing  length  turbulence 
models.  A  key  feature  of  the  present  calculations  is  that  grid  distributions 
which  provide  resolution  of  the  viscous  sublayer  regions  were  used,  as  are 
believed  important  for  three-dimensional  flows  of  this  type.  There  is  good 
overall  qualitative  agreement  between  the  computed  and  measured  flow 
structures.  The  location  of  the  initial  saddle-point  flow  separation  agrees 
well  with  the  experimental  flow  visualization.  The  strength  and  scale  of  the 
mean  flow  velocities  agrees  with  mt asurements ,  especially  within  the  reversed 
flow  region  upstream  of  the  leading  edge.  Computer-generated  particle  lines 
show  a  saddle-point  separation  followed  by  a  strong  downflow  toward  the 
endwall,  and  a  reversed  flow  region,  all  of  which  are  similar  to  numerous  flow 
visualization  experiments  for  horseshoe  vortex  flows.  The  computed  corner  flow 
has  a  multiple  vortex  structure.  There  are  quantitative  differences  between 
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computed  and  measured  velocities  in  the  corner  flow  region,  where  the  secondary 
vortex  Is  relatively  weak.  These  differences  are  likely  due  to  the  turbulence 
model  used,  and  possibly  numerical  error.  Computationally,  It  Is  significant 
that  these  solutions,  with  highly  stretched  meshes  to  provide  viscous  sublaver 
resolution,  required  only  100-200  iterations  for  convergence  using  a  60x50x40 
grid  (120,000  points). 
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